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A discussion of application of the Monte Carlo method to rarefied gas 
heat transfer is given. A sample problem of heat transfer through a rare- 
fied gas between infinite flat plates is treated. Hard sphere molecular col- 
lisions and a vail accomodation coefficient of 1 are assumed. The target 
molecule distribution is assumed to be Maxwellian. The results are com- 
pared to approximate analytical methods and to another Monte Carlo solution. 


IHTRODUCTION 

In rarefied gas transport problems the usual solutions using Navier 
Stokes* equation for momentum transfer or the Fourier equation for conduc- 
tion are no longer applicable. This is because these equations assume 
local isotropy and small gradients compared to the path lengths of the mole- 
cules. These assumptions are not valid in the case of rarefied gases. We 
must then resort to the more fundamental Boltzmann equation. This equation 
is difficult to treat by the usual analytical procedures because of its 
complexity. The Monte Carlo method allows us to reduce the complexity of 
the analysis at the expense of added numerical computation and is not re- 
stricted by the many simplifying assumptions generally made to allow analyt- 
ical solutions. 

The Monte Carlo procedure is a model sampling technique. We create a 
model and then follow histories of sample molecules through this model. 
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The sample histories are obtained by making choices at points of decision 
from the appropriate probability distribution. By averaging certain prop- 
erties of the sample molecules at various positions we can obtain the macro- 
scopic quantities of interest. 

For purposes of illustrating the Monte Carlo method let us consider 
the problems of heat transfer between parallel plates enclosing a rarefied 
gas. A discussion of much of the previous work on this problem is given in 
references 1 to 5. 

The model is shown in figure 1. The hot wall is considered to be at 
temperature T w q and the cold wall at temperature T w ^. A rarefied 
hard sphere molecule gas is contained between the walls. 

The sample molecule histories are started at the wall 0 by picking 
velocity components for the sample molecule leaving the wall from the ap- 
propriate distribution of velocities of the absorbed and re-emitted mole* 
cules. The space between the walls is divided into zones as shown in fig- 
ure 1. The sample molecule after leaving the wall will either pass through 
the first zone or have a collision with a target molecule in this zone. 

This will depend on whether the path length to collision for the sample 
molecule is longer or shorter than the distance the sample molecule must 
travel to pass through the zone. If the sample molecule passes through the 
first zone with no collision, it is started again at the beginning of the 
next zone with its velocity components unchanged. 

If there is a collision in the zone the point of collision is found at 
the end of the path length to collision. A target molecule is picked from 
the distribution of target molecule collision partners in the zone. The dis 
tribution of target molecule velocities is assumed uniformly Maxwellian 
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throughout each zone. A collision calculation is then carried out for the 
sample and target molecule collision partners, and new velocity components 
for the sample molecules are found. The history of the sample molecule is 
then continued from the point of collision with the new sample molecule 
velocity components. The target molecule collision partner is ignored after 
collision. 

A new path length to collision is found for the sample molecule, and 
this is compared in length to the distance from the point of collision to 
the next zone. If the path length to collision is greater, the sample mole- 
cule is started in the next zone with its velocity components unchanged. In 
the other case there is a collision in the zone and the collision calculations 
are repeated as before. 

If the sample molecule strikes the upper wall, it is assumed to be ab- 
sorbed and is reemitted with new velocity components picked from the ap- 
propriate distribution of velocities based on the upper wall temperature. 

After leaving the upper wall it is followed as before. The history is com- 
pleted when the sample molecule is incident on the lower wall and then a new 
sample history is begun. This is continued until the desired number of 
sample histories are completed. 

The transport and fluid characteristics of interest are the density 
distribution, temperature distribution, and heat transfer across the channel. 
These can be obtained as in reference 6 by locating scoring positions at 
various locations across the channel (see fig. l). By scoring the various 
characteristics of the sample molecules as they pass the scoring cross sec- 
tion, we can obtain the transport properties and fluid characteristics of 


interest. 
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The distribution of target molecules/ which is assumed uniform in each 
zone, will greatly affect the results since this will determine the path 
length to collision for the sample molecule and also the velocity component 
distribution of the target molecule collision partner in each zone. In the 
present analysis it was assumed that the target molecules in each zone were 
in a Maxwellian velocity distribution based on a different local temperature 
and density in each zone. Then by carrying out the Monte Carlo solution 
based on this assumed distribution of target molecules the temperature and 
density in each zone can be found. The Monte Carlo solution is then rerun 
using the new found local temperature and density in the Maxwellian distri- 
bution of target molecules. The problem is iterated in this manner until 
the density and temperature distribution found from the sample molecule 
histories agreed with the density and temperature distribution used for the 
target molecules. 

The assumption of a local Maxwellian distribution for the target mole- 
cules would be most applicable near equilibrium conditions, that is, cases 
with small temperature gradients and short mean free paths. It appears 
possible to extend the analysis to more realistic distributions of target 
molecules than were used in the present case as for instance a two-sided 
Maxwellian as used in an analytical solution in reference 5. 

In references 1 and 2 a Monte Carlo solution is used to treat a similar 
problem. In that analysis the solution is carried out with an assumed dis- 
tribution of target molecules in a tabular form. By scoring the velocities 
of many sample molecules as they pass through each zone they can obtain the 
velocity distribution of the molecules in each zone in a tabular form. 
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These are then used as the target molecule distributions in the next itera- 
tion. This process is continued until the sample molecule distribution 
found agrees with the target molecule distribution assumed. Then from the 
distributions in each zone the macroscopic quantities of interest can be 
found by numerical integration of the moments of the distribution. The pre- 
sent method, however, avoids the finding and storing of entire distributions 
of target molecule for each increment. 

In the Monte Carlo method since the distribution of target molecules is 
an approximation to the true distribution only conservation of mass is 
satisfied exactly since molecules are not allowed to "disappear" in transit 
between the walls. Conservation of momentum and energy can only be con- 
sidered in an average sense because the collision partners are not recorded. 

NOMENCLATURE 

C average thermal velocity (2KT/M) 1 / 2 

D channel height 

r x -u2 

erf(x) error function, (2/V*) Jq e du 

f probability distribution function 

f + ,f_ probability distribution function of molecules moving in positive, 

negative direction 

Kn Knudson number, Ty/D = M /-Jz Sp A D 

k Boltzmann constant 

M mass of molecule 

N number of sample molecules emitted from surface 0 in Monte Carlo 

run, proportional to the flux of molecules leaving surface 0 
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s + ,s_ 


pressure 

zone or scoring position number 
last scoring position number 
property of sample molecule 
averaged quantity Qf d^V 
averaged quantity f Qf + d 3 V 

random number between 0 and 1 
mutual collision cross section, na 2 

number of sample molecules through scoring position in positive, 
negative xg direction 
absolute temperature 
molecular velocity 
velocity after collision 


X]_,Xg,X 3 coordinates 


defined by eqs. (B6) and (B7) 
defined by eqs. (B6) and (B7) 

collision rate of sample molecule with target molecules 

path length to collision 

mean free path length 

dimensionless velocity, V/C 

defined by eq. (A5) 

mass density 

diameter of hard sphere molecule 

angle between sample molecule and target molecule velocities 
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Subscripts : 

A averaged 

p increment number 

p f last increment 

R relative velocity 

s sample molecule 

t target molecule 

w,0;w,l at wall 0, 1 

0,1 evaluated next to wall 0, 1 

1,2,3 coordinate directions 

+,- positive or negative direction 

ANALYSIS 


Start of Sample Molecule History at the Wall 
The sample molecule history will begin at the surface with temperature 
T w ^ 0 . To pick the components of velocity of the molecules leaving the sur- 
face the simplest assumption is that the molecules incident on the wall are 
perfectly accommodated, that is, are in a Maxwellian distribution at the 
wall based on the wall temperature. IMs assumption is discussed in refer- 
ence 3. The normalized Maxwellian velocity distribution of the molecules 
moving away from the wall (V 2 > 0) is given as 

Tf 2 f + ,0 d5v = ^3 (- |)] dVl dV 2 dv 3 V 2 >0 <*> 


The 

and 



/M is the number density of molecules moving away from the wall, 
is the thermal velocity T v,o) 1 ^ 2, If we consider the positive 
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Xg direction normal to the wall, the distribution of velocities of the mole- 
cules leaving the wall per unit time per unit area is Vgf + ^ 0 (refs. 1, 2, 
and 6). This can be transformed to cylindrical coordinates Vg = Vg, 

V x = V r cos 9, V 3 = V r sin 9 and normalized by (Vg) + Q = C q/V^ 

(ref. 1) to give 


f(e,v 2 ,v r ) = 


v 2 f +i0 d 3 v 

( V + , 0 


2V, 


r 


rtC 


0 


V £ 

exp I- -5 


V„ d9 dV P dV 
r 2 r 


f 9 f Vg f V r dG dV 2 dV r 


( 2 ) 


The distributions for fgfy fy 

2 r 


can be written separately as 


f 9 



(3) 



The velocity components of our sample molecule leaving the surface must be 
picked from these distributions. This same result was used in references 1 
and 6. A convenient way of picking from a distribution for the high speed 
computer is to transform the distribution to a uniform distribution in R 
by setting the random number R equal to the cumulative distribution func- 
tion. For instance 



( 4 ) 


Then by using a high speed computer to generate a random number R 
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between 0 and 1, we can obtain 0 from equation (4) such that for a 
large number of samples picked in this manner the distribution in equa- 
tion (3a) will be satisfied. 


Similarly we can pick Vg from 



or, since picking R is equivalent to picking 1 - R, 

Vg = (-Cq In (6) 

so that by picking Ry we can obtain Vg from equation (6.). The velocity 

U 

V r is obtained in the same manner and is given by 

V r = (-C§ In Ry^ 1 / 2 (7) 

The Vg, V r , and 0 picked then give the direction and velocity of the 
sample molecule as it leaves the wall. This same result was used in refer- 
ences 1 and 6. 

The sample molecule is then followed through the first zone until 
either a collision occurs with another molecule or the sample molecule 
passes through the zone. This will depend on whether the path length of the 
sample molecule to collision is shorter or longer than the path length 
through the zone. 

Calculation of the Path Length A to Collision 
for the Sample Molecule 

The probability that a sample molecule will collide in the incremental 
path length A to A + dA is given in reference 7 as 


10 


exp 


(-V^,s) 


( 8 ) 


|i,S 


where A,, _ is the mean free path to collision of the sample molecule mov- 
s 

ing at velocity V s in that zone. 

We can then pick a path length to collision for the sample molecule from 
this distribution by using the same procedure described earlier 


^ - “\i,s ■*' n 


(9) 


V,s 


To use this relationship, we must first calculate the mean free path ^ 
in the zone. 

As shown in appendix A the mean free path of a sample molecule moving 
at velocity V s through a Maxwellian gas at density p with a thermal 
velocity C is (eq. (A8)) 

2 1 / 2 Kn u 


|i,S 


_P_ 

pa 


77T- + <erf M » ) (“■ + at) 


( 10 ) 


where p s is the velocity of the sample molecule nondimensionalized by the 
thermal velocity of the zone (p g = V s /C), D is the distance across the 
channel and Kn is the usual definition of Knudson number for hard sphere 
molecules in a Maxwellian distribution 



M 

V2 Sp A D 


( 11 ) 


Equation (10) can then be used in equation (9) to obtain the distance to col 
lision for the sample molecule in the zone. If there is a collision in the 
zone, the next step is to pick a collision partner and calculate the new 
direction and velocity of the sample molecule after collision. 
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Calculating New Sample Molecule Velocity Components 
and Direction After Collision 

The distribution of target molecules hit by the sample molecule is 
given in appendix B following the derivation in reference 7. From this dis- 
tribution we pick a target molecule collision partner. The new velocity 
of the sample molecule after collision with the target molecule is given 
following the derivation in reference 1. 

V £ s - \ ( v 2s + Vgt) + \ V R (1 - 2* 2 ) 

V L - | ( V l 8 + V lt> + T r( 1 * b2 ) l/2 H 

V ! s * 2 t V 3s + V 3t) + V R (1 - t,2) l/2 a 

where 

V R = [(Vis " V i t ) 2 + (V 2s “ V 2t ) 2 + (V 3S - V 3t ) 2 ] l/2 

The values of H and H are obtained as discussed in reference 1 by pick- 
ing two rapdom numbers and using them in the following equations 
H = 2Rh -: 1, H = 2R- - 1, and b 2 = H 2 + H 2 (13) 

where b 2 must be less than one. If b 2 is greater than one, a new set of 
random numbers must be chosen to find H and H. 

Completion of Sample Molecule Histories 
After collision a new path length to collision is found and this is 
compared to the distance the sample molecule must travel in its new direc- 
tion to leave the zone. If this distance is smaller than the new path length 
to collision, the molecule is then started through the next zone as shown in 
figure 1. This process is continued until the sample molecule returns to the 




( 12 ) 
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wall 0 and then a new sample molecule history is started. This process 
is continued for the desired number of sample histories. 

Scoring to Find the Macroscopic Flow Properties 
The macroscopic fluid characteristics we must obtain are the density 
distribution and temperature distribution across the channel since these are 
needed in the target molecule distribution. Also of interest is the net 
heat transferred across the channel. These properties are obtained as fol- 
lows. Scoring positions are located at various distances across the channel 
as shown in figure 1. The average quantity of Q transported across the 
scoring cross section p in the positive Xg direction by the sample mole- 
cules crossing it can be written as 



where S + ^ is the number of sample molecules passing across the scoring 
cross section p in the positive direction and Q is some quantity 

each sample molecule carries. Similarly the average quantity Q trans- 
ported in the negative X2 direction is 



(15) 


Since there is no net flow across the channel and since all sample histories 
start and end at the wall 0., 


S 


+>.P 



(16a) 


and 
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If V+ 



We can then write 


p(V2Q) 


+,P 


£_ 

+,P 



(16b) 


(17) 


The number of sample molecules that pass the scoring cross section at p 
in the positive V 2 direction S + ^ divided by the toted, number of sample 
molecule histories started at wall 0, N, can be related to the mass flux 
passing in the positive V 2 direction at p by 


ti2 


N 


< 0 ( <0 


where p, 0 (v), n the mass flux leaving wall 0 is equal to 

• y U * y U 

(ref* 1), Combining equations (17) and (18) results in 


p(V 2 Q) p 



(18) 

p+,o c +,oA l/2 


(19) 


If Q is taken as l/v 2 , equation (19 ) becomes 



The average density in the channel is then obtained by averaging the density 


of all the scoring cross sections from p = 0 to p = p^. 
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P +,0 P f 


P 0 M Pf 


V 1 

2 

P=1 


+,0 


As shown in reference 1 the collision solution (Kn -*• «°) is given by 

\l/2 

V 

~ 'Kn-** 0 


f— ^ 

\ P +,0/ f 


+ 1 


\^w , 0/ 

/ T .\ 1/2 

I 1 

yv,o) 


Similarly we can obtain the local pressure from 


(21a) 


(21b) 


1 <1 + V r> 


< pv 2 > 


3 P / «ir2 \ _ ~^r -Wfl-1/2 


+ ^ Q (pV /+ ;0 3C + ^ 0 Nrt' 

P. ,C 2 



i + v l 


i 


i + v 2 


( 22 ) 


since (pVg^ q = + ;Q + ;Q as given in reference 1 and V 2 = + v |. 

Dividing equation (22) by p/p + as obtained from equation (20), then gives 
the local temperature 


T <V 1 + T ? > < P V 2> + ,0 


(23) 


3C 


± 2 ° 


J +,0 


as shown in reference 1 the collisionless solution (Kn -*• °°) is given by 


(j-\ . fei 

y^oyKn-* 00 V-w, 


’••r 


( 24 ) 


Finally, the heat transfer across the channel can be obtained from 
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p <(vl + vg)) . P+; i/2 ;9 

^ EjA^N 


^+,p 

X 


-,P 


(v| + V^) - 



(Vi + V2) 


(25) 


This can be nondimens ionali zed by = — (p(Vg(v| + v|))) + Q = - A j (p + Q C^ Q ). 

the heat leaving the vail as given in (ref. l) and then divided by P + /p^ 
to give 


A/2 


^ 0. - p 

P A C^0 21,0 


+,0 


+ >p ->p 

(vf + vf) - (vl + v|) 


(26) 


This can be compared to the collisionless result given in reference 1 



In this manner the local temperature and density in each zone can be obtained. 
This is compared to the assumed temperature and density distribution used in 
the local Maxwellian distribution of the target molecules in each zone and 
the results are iterated till agreement is obtained. 


RESULTS AND CONCLUSIONS 


The results for the density distributions and temperature distributions 
are shown in figures 2 and 3. The heat- transfer results are shown in 
figure 4. These results were compared with the Gross-Ziering eight moment 
half range solution as given in reference 4 and Haviland and Lavin nonlinear 
four moment solution given in reference 5 and their Monte Carlo solution 
given in reference 2. 

Both Monte Carlo solutions are in better agreement with the Gross- 
Ziering solution and poor agreement with the nonlinear four moment solution. 
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The present technique is very flexible and can be readily extended to 
different types of molecular collisions and to different boundary condi- 
tions. Using the Maxwellian distribution for the target molecules resulted 
in a nonconstant heat transfer across the channel. This was averaged to 
obtain the present result. More realistic target molecule distributions , 
however, can be used in the analysis. An important shortcoming of the 
present method is the need of a high speed computor and extensive comput- 
ing time. Each solution required about 1 hour of running time on the IBM 
7094. 
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APPENDIX A 

EVALUATION OF MEAN FREE PATH 


The number of collisions per unit time d0(Vjj,Vj.) of a sample molecule 
moving at velocity V g through target molecules in velocity volume space 
d^Vj. for hard sphere molecules is given in reference 6 as 
d0 = pf t V R S d 3 ^ (Al) 

•where V R is the velocity of the target molecules relative to the sample 
molecule velocity before collision, V R .^ = ^t,i “ ^s,!’ as s ^ own i- n fig* 

p 

ure 5, and S is the mutual collision cross section Jto where a is 

the diameter of the molecule. The relative velocity V R as seen from fig- 
ure 5 can be written as 

V R = (vf + V§ - 2V s V t cos cp* ) 1 / 2 (A2) 

where cp* is the angle between and V g . For the case of a Maxwellian 

distribution following the analysis given in reference 7 the collision rate 
can be written as 


de(p t ,qp',0) = 


pps 

^75 



sin cp* dcp' d0’ 


d^ 


(A3) 


where p is the nondimensionalized velocity V/C and the angles 9' 
and 0' are shown in figure 5. To obtain the total collision rate for the 
target molecules over all velocities, this must be integrated over 9', 0', 
and |i, . We can integrate this over 0' from 0 to 2 k and integrate 
over 9' from 0 to n to give 


ae(nt) = ^J72 [ ex P (-4)] Pt u R,A 


(A4) 


where 
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or 


,a = f + m§ - 2^ t cos q') 1 / 2 S 1 g ~ 

Jq 


dcp' 


R,A 


\ + s if H 

I Wt + ^s/ 3 ^t if > 


(A5) 


We can then integrate over | from 0 to 00 to obtain the total col- 
lision frequency for a sample molecule moving at velocity p s through a 
Maxwellian gas as 


0 = 


pSC 

[exp (-Pg) 

M 

( «V2 


+ (erf |i s ) 


(^)j 


(A6) 


The dimensionless mean free path is then found as discussed in reference 7 
by dividing the dample molecule velocity by the total collision rate to 
give 

d s M 


^S_Vg 

D ~ 0 


DpS ^^T 72 [ 6XP (_M s)) + < erf ^b) Y s + 2^)} 


(A7) 


Using the definition of Knudson number as given by equation (11) enables us 
to write 




= Kn2 


. 1 / 2 , 


1 j" (exp (-n|j) 

J72 


<«fP=; 


+ (erf n_) 


f s + 2O. 


-1 


(A8) 
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APPENDIX B 

PICKING THE TARGET MOLECULE COLLISION PARTNER 


The velocity distribution of target molecules that the sample molecule 
will collide with is discussed in reference 1 is given by 
deC^jCp'jd) [exp (-p|)]ii|p R sin cp’ dq>’ d 0’ d^ 

6 ntexp (-p|)3+ * 3 / 2 (erf p s )[x s + j 

The distribution in 0' for the target molecules is readily seen to be 


f e , d0' 


dQ 1 

2it 


(B2) 


Then the 0' can be picked from this distribution by 

9' = 2*R e (B3) 

The distribution of for the target molecules is obtained from the 

marginal distribution 


4[exp (-n|) J n||ip A dp t 
^ * [[exp (-nj)] + « l/2 (erf d g )(d s + 

We can then pick from this distribution as before 

r(n s ^) + 

t [exp (-M-g ) 1 + Jt^/ 2 (erf fT7[j” + 
where when |a_ > 

8 ^ 3 

r(d s ,d t ) = - 2 u s d t |^exp(-^)] + d s V« ( erf *0 " I vj [sxp(-u|)J 

- ^[ ex P (-P?)] + (erf 

T1 = 0 


(B4) 


(B5) 


(B6) 


and when u < u + 

S u 
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Then components of velocity of the target molecule must still he 

transposed to the coordinate system of the channel. Taking in the 

same direction as the coordinate as shown in figure 5 and V^. in the 

same direction as the coordinate we can obtain the component of 

velocity for the target molecule in the primed system as 

(a 1 = p , cos cp' 

2t t 

(i^ = (i^. sin cp' cos 9 ' 

l-i 1 = pi, sin cp' sin 9 ' 

3t t 

Then by a simple rotation of the coordinate system through an angle cp 
around the V3 axis the components of velocity of the target molecule in 
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the channel coordinates are given by 
^3t - ^3t 

U 2t = ^ cos «p - sin cp 

m . = n’ sin cp + u* cos cp 
-i-t 2t It 


> 


(B12) 
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Figure 3. - Temperature distribution across channel Knudson number, 2. 



Figure 4. - Heat transfer across channel. 



Figure 5. - Velocity of sample and target molecules at collision. 





